/*

class  implementation

*/


#include <TerrainGuarding.h>

extern bool PARALLEL;
extern bool USE_THRDS;

TGInstance::TGInstance(std::string& input_file, input_mode input_type) : input_type(input_type)
{
	
	index i,j;
	size_t m,n;

		
	std::ifstream infile;
        infile.open(input_file.c_str());


        if (!infile.is_open()) { std::cout << "error: error reading input file." << std::endl; std::exit(1);}

	

	// input problem via ILP

	
	infile >> m >> n >> D;

	_num_points = m; _num_guards = n;

	// locate data to store visibilty matrix
        // memory alocation for arrays 

   	
     V.allocate(m,n,0);
	 demands.allocate(m,0);
	 weights.allocate(n,0);

	long readings = 0;
	bool val;

	while(++readings<=D)
	{
		infile >> i >> j >> val;
		V[i][j] = val;			
	}

	readings = 0; real x; 

	// check instance
	inst_type t1 = GUARDING,t2 = GUARDING;
	while(++readings<=m) {infile >> x; demands[readings-1] = x; if(x>1.0) t1=MULTI_GUARD;  }; 
	readings = 0;
	while(++readings<=n) {infile >> x; weights[readings-1] = x; if(x>1.0) t2=WEIGHT_GUARD;  };


	infile.close();


	if(t1 == MULTI_GUARD && t2==GUARDING) instance_type = MULTI_GUARD;
	if(t1 == GUARDING && t2==WEIGHT_GUARD) instance_type = WEIGHT_GUARD;
	if(t1 == GUARDING && t2==GUARDING) instance_type = GUARDING;
	if(t1 == MULTI_GUARD && t2==WEIGHT_GUARD) instance_type = WEIGHT_MULTI_GUARD;


	
	
		
	  
}



input_mode TGInstance::getInputType() const { 	return this->input_type; }

inst_type TGInstance::getInstType() const {return instance_type;}


std::ostream& operator << (std::ostream& out, const TGInstance& T)
{
	  // check readings
	out << "guards: " << T._num_guards << " points: " << T._num_points << " non-zeros:" << T.D << std::endl;
        out << T.V << std::endl;
	out << T.weights << std::endl; 
	out << T.demands << std::endl;

	return out;

}

bool& TGInstance::update_visible(index p, index q) { return V[p][q]; }

bool TGInstance::visible(index p, index q) const { return V[p][q]; }

bool TGInstance::left_visible(index p, index q) const { return ( (p<q) && visible(p,q)); }

bool TGInstance::right_visible(index p, index q) const { return ( (p>q) && visible(p,q)); }



size_t TGInstance::num_guards() const { return _num_guards;}
size_t TGInstance::num_points() const { return _num_points;}
real TGInstance::weight(index g) { return weights[g];}
real TGInstance::demand(index p) { return demands[p];}


/*
	TGSolver class implementation

*/


TGSolver::TGSolver(TGInstance& T, double omega) : T(T), omega(omega), cost(0), d_min(1)
{
	
	std::cout << " ******************* 1.5D Terrain Guarding Solver ****************** " << std::endl;
	
	Timer compute_time;
	
	compute_time.start();
	
	status_code status;
	
	if(T.getInputType() == GEOMETRIC) std::cout << "Sorry, not yet implemented." << std::endl;
	if (T.getInputType() == ILP) std::cout << " -> using ILP formulation for 1.5D Terrain Guarding problem as input:" << std::endl;
	std::cout << "\r-> points: " << T.num_points() << " guards: " << T.num_guards() << 
		" instance type: " << ((T.getInstType()==WEIGHT_GUARD || T.getInstType()==GUARDING ? "weighted-version" : "version with demands")) << std::endl;
	std::cout << "\r-> error: " << omega << std::endl;
	std::cout << "---------------------------------------------------------------------" << std::endl;	
	status = LPSolver(omega);
	
       
	status = RoundGuards(1/5.0);
	
	if(status!=ROUND_FIN) 
	{
	   bool useThreads = (tthread::thread::hardware_concurrency()>=2) && USE_THRDS;
	  if(useThreads) 
	  {
		std::cout << "\r-> using 2 threads to solve left and right guarding problem in parallel." << std::endl;
		tthread::fast_mutex mutex;
		Arg argL(T,left_demand,guards,cost,mutex);
		Arg argR(T,right_demand,guards,cost,mutex);
		tthread::thread left_guard_thread(LeftGuardingThreadCode,(void*) &argL);
		tthread::thread right_guard_thread(RightGuardingThreadCode, (void*) &argR);
		
		left_guard_thread.join();
		right_guard_thread.join();
		
		status = FIN;
		
	  }
	  else{
	  status = LeftGuarding();
	  status = RightGuarding();

	 
      }
     }
	
	
	if(status==FIN || status==ROUND_FIN) 
	{
		compute_time.stop();

		std::cout << "-----------------------------------------------------------------------" << std::endl;
		std::cout << "*** guard cost:   " << total_cost() << std::endl;
		std::cout << "*** compute time: " << compute_time.secs() << std::endl;
		std::cout << "-----------------------------------------------------------------------" << std::endl;
		std::cout << std::endl << std::endl;
		}
	else
	{	
		std::cout << "-----------------------------------------------------------------------" << std::endl;
		std::cout << "Error in computing the solution! " << std::endl;
		std::cout << std::endl << std::endl;
	}
}

status_code TGSolver::LPSolver(dreal omega)
{
    LP_status status;
    if(omega>0)
	{
		if(VERBOSE) std::cout << "solving left guarding problem ... " << std::endl;
		// use packing formulation
		size_t m = T.num_points(); size_t n = T.num_guards();		

		bool  *A = T.V.CArrayPtr();  // A^T in column-wise storage
		real *b = T.weights.CArrayPtr();
		real *c = T.demands.CArrayPtr();

                
		// packing solution
		vector<dreal> pack_sol(m,0);

		// covering solution (needed for the algorithm)
		cover_sol.allocate(n,0);

		// solved as fractional packing-covering problem

		
		bool parallel_treshold =(T.num_points() > 400);
		
		#ifdef __CUDACC__
			if(parallel_treshold && PARALLEL) 
			{
				std::cout << "-> method: GPU based approximate packing-covering solver " << std::endl; 
				status = cuGKPCSol(A,b,c,pack_sol.CArrayPtr(),cover_sol.CArrayPtr(),n,m,omega); 
			}
			else 
			{
				std::cout << "-> method: CPU based approximate packing-covering solver" << std::endl;
				status = GKPCSol(A,b,c,pack_sol.CArrayPtr(),cover_sol.CArrayPtr(),n,m,omega);  
			}
		#else 
			status = GKPCSol(A,b,c,pack_sol.CArrayPtr(),cover_sol.CArrayPtr(),n,m,omega); 
		#endif	

		
       }

	if(status==LP_FEAS) return LP_SOL;
	else return LP_NON_SOL;
}


inline void TGSolver::decompose(index p,vector<dreal>& x, dreal& left_frac, dreal& right_frac) 
{
	integer n = T.num_guards();
	left_frac = right_frac = 0.0;

	if(T.demand(p)>0)  
	{  // if p not satisfied in rounding, define left and right demand
	
		for(index g=0;g<n;g++)
		{ 
			left_frac  += x[g]*T.left_visible(g,p);
			right_frac += x[g]*T.right_visible(g,p);
		}

	
		left_frac = (1+1.0/d_min)*(left_frac+0.5*x[p]); right_frac = (1+1.0/d_min)*(right_frac+0.5*x[p]);
	
		if(T.getInstType()==WEIGHT_GUARD || T.getInstType()==GUARDING)
		{
			left_frac = (left_frac>1) ? 1 : 0;
			right_frac = (right_frac>1) ? ( (left_frac<1) ? 1 : 0 ) : 0;
		}
		else 
		{
			left_frac = ceil(left_frac);
			right_frac = ceil(right_frac);
		}
      }
      else {left_frac=0; right_frac=0; }
		
}

status_code TGSolver::RoundGuards(dreal alpha)
{
	std::cout << "\r-> rounding large fractional value guards ..." << std::endl;
	
	// initialize guards
	guards.allocate(T.num_guards(),0);
	left_demand.allocate(T.num_points(),0);
	right_demand.allocate(T.num_guards(),0);	
	cost = 0.0;

	//if(VERBOSE) {std::cout << "fractional solution :" << cover_sol << std::endl;}
	
	// round point guards
	if(VERBOSE) std::cout << "rounded guards: " << std::endl;
	for(index g=0;g<T.num_guards();g++)
	{ 
		if(cover_sol[g]>=alpha)
		{ 
			
			if(VERBOSE) std::cout << g << " ";

			// remove large frac. value guards, update solution
			guards[g] = 1; cost += T.weight(g); cover_sol[g] = 0.0;

			for(index p=0;p<T.num_points();p++)
			{ 							
				T.demands[p] = ( (T.demands[p]>0) && T.visible(g,p)) ? (T.demands[p]-1) : T.demands[p];
				// T.update_visible(p,g) = T.update_visible(g,p) = 0;
			}
			
		}
				
	
	}

	if(VERBOSE) std::cout << "\n\rcost after rounding: " << cost << std::endl;

		
	// find minimum demand
	d_min = std::numeric_limits<int>::max();
	integer d_max = std::numeric_limits<int>::min();

	for(index p=0;p<T.num_points();p++) 
	{
		d_min = (T.demand(p)>0) ? min<integer>(d_min,T.demand(p)) : d_min; d_max = max<integer>(d_max,T.demand(p)); 
		
	}  
	
	       
	 if(d_max==0) return ROUND_FIN; // check if guard rounding finishes the algorithm

	// define left and right guarding problem
	
	
	for(index p=0;p<T.num_points();p++)
	{
	  dreal d_L,d_R;
	  decompose(p,cover_sol,d_L,d_R);
	  
	  left_demand[p] = d_L;
 	  right_demand[p] = d_R;
	  
	 	  
    }  
	

	if(VERBOSE){
		// std::cout << "\nupdated instance: " << std::endl;
		// std::cout << T << std::endl;
		// std::cout << "\nleft, right residual demand:" << std::endl;
		// std::cout << left_demand << std::endl;
		// std::cout << right_demand << std::endl;
	}

	return FIN;
};


status_code TGSolver::LeftGuarding()
{
		
	integer m = T.num_points();
	integer n = T.num_guards();
	
	
	
	std::cout << "\r-> solving left guarding problem ...";

	if(T.getInstType() == WEIGHT_GUARD || T.getInstType() == GUARDING)
	{
		std::cout << "weighted version" << std::endl;
		vector<integer> trigger_guard(m,NOT_SET);
		vector<integer> cnt_left_guards(m,0);
		vector<real> w(m,0); // residual weight

		for(integer g=0;g<n;g++) w[g] = T.weight(g); // copying weights

			
		real w_min = INFTY;
		integer g_p = m;
	

		// finding left guards for any point p
		for(integer p=0;p<m;p++)
		{
			if(cnt_left_guards[p] > 0) continue;  // consider only not covered points and are to be guarded from left
			if(left_demand[p]>0)
			{ 
				for(integer g=0;g<p;g++) if(T.visible(g,p) && (w_min>w[g])) {w_min = w[g]; g_p = g; } // find min-weight-guard
				// guard for point p is g_p 
				trigger_guard[p] = g_p; 
				cost += T.weight(g_p); guards[g_p]=1;
				// update residual weights
				for(integer g=0;g<p;g++)  w[g] -= w_min*T.visible(g,p);
				// which points are also covered with g_p ?
				for(integer q=g_p+1;q<m;q++) cnt_left_guards[q] += T.visible(g_p,q);
				
		
			}
			
			
		}
		

		// pruning step: remove redundant guards
		for(integer p=m-1;p>=0;p--)
		{			
			if(trigger_guard[p]==NOT_SET) continue;
			if(cnt_left_guards[p]>1)
			{ 
				g_p = trigger_guard[p];
				for(integer q=g_p+1;q<m;q++) cnt_left_guards[q]-= T.visible(g_p,q);
				trigger_guard[p] = NOT_SET; 
				cost -= T.weight(g_p); guards[g_p] = 0;
			}
		}

				
	}
	


	// not implemented yet

	else if (T.getInstType() == MULTI_GUARD)
	{
		integer m = T.num_points();
		integer n = T.num_guards();
		
		std::cout << "demand version" << std::endl;
		vector<integer> trigger_guard(m,NOT_SET);
	
	
	
		std::cout << "\r-> solving left multi-guarding problem ...";
		
		for(integer g=0;g<n;g++)
		{
			for(integer p=m-1;p>g;p--) if(T.visible(g,p) && left_demand[p]>0){ left_demand[p]--; guards[g]=1;}
		}
	
		
	
	}

	return FIN;
}


status_code TGSolver::RightGuarding()
{
	// symmetric version of LeftGuarding algorithm
	integer m = T.num_points();
	integer n = T.num_guards();

	std::cout << "\r-> solving right guarding problem ...";

	if(T.getInstType() == WEIGHT_GUARD || T.getInstType() == GUARDING)
	{
		std::cout << "weighted version" << std::endl;
		vector<integer> trigger_guard(m,NOT_SET);
		vector<integer> cnt_right_guards(m,0);
		vector<real> w(m,0); // residual weight

		for(integer g=0;g<n;g++) w[g] = T.weight(g); // copying weights

			
		real w_min = INFTY;
		integer g_p = m;

		

		// finding right guards for any point p
		for(integer p=m-1;p>=0;p--)
		{
			if(cnt_right_guards[p] > 0) continue;  // consider only not covered points and are to be guarded from right
			if(right_demand[p]>0)  // should be covered from right?
			{ 
				for(integer g=m-1;g>p;g--) if(T.visible(g,p) && (w_min>w[g])) {w_min = w[g]; g_p = g; } // find min-weight-guard
				// guard for point p is g_p 
				trigger_guard[p] = g_p;
				cost += T.weight(g_p); guards[g_p]=1;
				// update residual weights
				for(integer g=m-1;g>p;g--)  w[g] -= w_min*T.visible(g,p);
				// which points are also covered with g_p ?
				for(integer q=g_p-1;q>=0;q--) cnt_right_guards[q] += 1*T.visible(g_p,q);
			}
			
		}
		
		
		// pruning step: remove redundant guards
		for(integer p=0;p<m;p++)
		{			
			if(trigger_guard[p]==NOT_SET) continue;
			if(cnt_right_guards[p]>1)
			{ 
				
				g_p = trigger_guard[p];	
				for(integer q=g_p-1;q>=0;q--) cnt_right_guards[q]-=1*T.visible(g_p,q);
				trigger_guard[p] = NOT_SET; 
				cost -= T.weight(g_p); guards[g_p] = 0;				
			}
		}

		
		
	}
	


	// not implemented yet

	else if (T.getInstType() == MULTI_GUARD)
	{
		integer m = T.num_points();
		integer n = T.num_guards();
		
		std::cout << "demand version" << std::endl;
		vector<integer> trigger_guard(m,NOT_SET);
	
	
	
		std::cout << "\r-> solving left multi-guarding problem ...";
		
		for(integer g=m-1;g>=0;g--)
		{
			for(integer p=0;p<g;p++) if(T.visible(g,p) && right_demand[p]>0){ right_demand[p]--; guards[g]=1;}
		}
	
		
	
	}

	return FIN;
	
}


dreal TGSolver::total_cost() const {return cost;}



/*

Thread enabled left guarding right guarding

*/

Arg::Arg(TGInstance& T,vector<index>& demand, vector<bool>& guards, dreal& cost, tthread::fast_mutex& mutex) : T(T), demand(demand), guards(guards), cost(cost), mutex(mutex) {}
		
	




void LeftGuardingThreadCode(void *args)
{
		Arg* arg = (Arg*) args;
		
		integer m = (arg->T).num_points();
	    	integer n = (arg->T).num_guards();
	
	    
		vector<integer> trigger_guard(m,NOT_SET);
		vector<integer> cnt_left_guards(m,0);
		vector<real> w(m,0); // residual weight

		for(integer g=0;g<n;g++) w[g] = (arg->T).weight(g); // copying weights

			
		real w_min = INFTY;
		integer g_p = m;
	

		// finding left guards for any point p
		for(integer p=0;p<m;p++)
		{
			if(cnt_left_guards[p] > 0) continue;  // consider only not covered points and are to be guarded from left
			if(arg->demand[p]>0)
			{ 
				for(integer g=0;g<p;g++) if((arg->T).visible(g,p) && (w_min>w[g])) {w_min = w[g]; g_p = g; } // find min-weight-guard
				// guard for point p is g_p 
				trigger_guard[p] = g_p;

				// use locking mechanism
				arg->mutex.lock();
				arg->cost += (arg->T).weight(g_p); arg->guards[g_p]=1;
				arg->mutex.unlock();
				
				// update residual weights
				for(integer g=0;g<p;g++)  w[g] -= w_min*(arg->T).visible(g,p);
				// which points are also covered with g_p ?
				for(integer q=g_p+1;q<m;q++) cnt_left_guards[q] += (arg->T).visible(g_p,q);			
		
			}		
			
		}	

		// pruning step: remove redundant guards
		for(integer p=m-1;p>=0;p--)
		{			
			if(trigger_guard[p]==NOT_SET) continue;
			if(cnt_left_guards[p]>1)
			{ 
				g_p = trigger_guard[p];
				for(integer q=g_p+1;q<m;q++) cnt_left_guards[q] -= (arg->T).visible(g_p,q);
				trigger_guard[p] = NOT_SET; 
				
				// enable lock
				arg->mutex.lock();
				arg->cost -= (arg->T).weight(g_p); arg->guards[g_p] = 0;
				arg->mutex.unlock();
				// lock disabled
			}

		}

				
}



void RightGuardingThreadCode(void *args)
{
		// symmetric version of LeftGuarding algorithm
		Arg* arg = (Arg*) args;
		
		integer m = (arg->T).num_points();
	    integer n = (arg->T).num_guards();

		
		vector<integer> trigger_guard(m,NOT_SET);
		vector<integer> cnt_right_guards(m,0);
		vector<real> w(m,0); // residual weight

		for(integer g=0;g<n;g++) w[g] = (arg->T).weight(g); // copying weights

			
		real w_min = INFTY;
		integer g_p = m;

	

		// finding right guards for any point p
		for(integer p=m-1;p>=0;p--)
		{
			if(cnt_right_guards[p] > 0) continue;  // consider only not covered points and are to be guarded from right
			if(arg->demand[p]>0)  // should be covered from right?
			{ 
				for(integer g=m-1;g>p;g--) if((arg->T).visible(g,p) && (w_min>w[g])) {w_min = w[g]; g_p = g; } // find min-weight-guard
				// guard for point p is g_p 
				trigger_guard[p] = g_p;
				
				// enable lock
				arg->mutex.lock();
				 arg->cost += (arg->T).weight(g_p); arg->guards[g_p]=1;
				arg->mutex.unlock();
				// update residual weights
				for(integer g=m-1;g>p;g--)  w[g] -= w_min*(arg->T).visible(g,p);
				// which points are also covered with g_p ?
				for(integer q=g_p-1;q>=0;q--) cnt_right_guards[q] += (arg->T).visible(g_p,q);
			}
			
		}
	
	
		// pruning step: remove redundant guards
		for(integer p=0;p<m;p++)
		{			
			if(trigger_guard[p]==NOT_SET) continue;
			if(cnt_right_guards[p]>1)
			{ 
				
				g_p = trigger_guard[p];	
				for(integer q=g_p-1;q>=0;q--) cnt_right_guards[q]-=1*(arg->T).visible(g_p,q);
				trigger_guard[p] = NOT_SET; 
				
				// enable lock
				arg->mutex.lock();
				arg->cost -= (arg->T).weight(g_p);  arg->guards[g_p] = 0;	
				arg->mutex.unlock();
				// disable lock
			}
		}

	
	
};








